#include <algorithm>
#include <iostream>
#include "2pdm_optimal_irreps.h"
#include "npdm_typedefs.h"

//----------------------------------------------------------------------------------------------------------------------------------------------------------
// Constructor for 2PDM_index object.

TwoPdmBlueprint::TwoPdmBlueprint (int kmax) {

  std::vector< std::vector< std::tuple<int,CRE_DES,ALPHA_BETA> > > ij_vec, lm_vec;
  kmax_ = kmax;

  // Build creation 2-indices ordered by site
  get2opsPerms(ij_vec,CREATION);
  // Build destruction 2-indices ordered by site
  get2opsPerms(lm_vec,DESTRUCTION);

  printAnalyticalResults();

  buildIrrepOpStrings( ij_vec, lm_vec );
  buildSpinDMRGmap();

}

//----------------------------------------------------------------------------------------------------------------------------------------------------------

void TwoPdmBlueprint::printAnalyticalResults()
{
   // Number of sites or spatial orbitals
   int n = kmax_;

   int n2 = n*n;
   int n3 = n2*n;
   int n4 = n3*n;

//   int total = 3*n4/4 - n3/2 + 5*n2/4 - n/2;
   int total = (n*(3*n3 - 2*n2 + 5*n -2))/4;
   std::cout << "Irreducible spin-strings = " << total << std::endl;

}

//----------------------------------------------------------------------------------------------------------------------------------------------------------
// Build ordered non-zero permutations of 3 creation/des ops.
// Sorted by site-index

void TwoPdmBlueprint::get2opsPerms(std::vector< std::vector< std::tuple<int,CRE_DES,ALPHA_BETA> > > & pq_vec, CRE_DES cre_des) {

  std::vector< std::tuple<int,CRE_DES,ALPHA_BETA> > ops(2, std::make_tuple(0,CREATION,ALPHA) );
  std::vector< std::tuple<int,CRE_DES,ALPHA_BETA> > tmp;

  ALPHA_BETA ab;
  for (int i=-kmax_; i<=kmax_ ; i++) {
    if (i==0) continue;
    if ( abs(i)/i == 1) ab = ALPHA;
    else ab = BETA;
    ops[0] = std::make_tuple(abs(i), cre_des, ab);
    
    for (int j=i+1; j<=kmax_ ; j++) {
      if (j==0) continue;
      if ( abs(j)/j == 1) ab = ALPHA;
      else ab = BETA;
      ops[1] = std::make_tuple(abs(j), cre_des, ab);
      tmp = ops;
      std::sort( tmp.begin(), tmp.end() );
      pq_vec.push_back( tmp );

    };
  };

}

//----------------------------------------------------------------------------------------------------------------------------------------------------------

